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The ordering of the classical Heisenberg model on the pyrochlore lattice with the antiferromag- 
netic nearest-neighbor interaction Ji and the ferromagnetic next-nearest-neighbour interaction J2 
is investigated by means of a mean-field analysis and a Monte Carlo simulation. For a moderate 
J2/Ji-value, the model exhibits a first-order transition into an incommensurate multiple-g ordered 
state where multiple Bragg peaks coexist in the spin structure factor. We show that there are two 
types of metastable multiple-g states, a cubic symmetric sextuple-g state and a non-cubic symmet- 
ric quadruple-q state. Based on a Monte Carlo simulation, we find that the cubic sextuple-g state 
appears just below the first-order transition temperature, while another transition from the cubic 
sextuple-5 state to the non-cubic quadruple-5 state occurs at a lower temperature. 

PACS numbers: 75.10.Hk, 75.30.Kz, 75.50.Ee, 75.40.Mg 



I. INTRODUCTION 

Recently, geometrically frustrated niagnets have at- 
tracted much interest due to their unconventional or- 
dering behaviors [iHSJ. Spin systems on the pyrochlore 
lattice, which consists of a three-dimensional network 
of corner-sharing tetrahedra (FiglT]), are typical exam- 
ples of such geometrically frustrated magnets. The 
classical Heisenberg magnet with the antiferromagnetic 
nearest-neighbor (NN) interaction is known to exhibit 
no magnetic long-range order even at zero temperature 
[3-0I ■ This is due to the macroscopic degeneracy of the 
ground state induced by geometrical frustration. Since 
such a high degeneracy is realized via a fine balance 
among frustrated interactions, it might be lifted by in- 
troducing small perturbations, e.g. the further- neighbor 
interactions |8l-[l3| . the quen ched randomness |13l - [l7| . or 
the lattice distortion [l8l - l20| . The lifting of the degener- 
acy as a result of small perturbations might lead to an 
exotic magnetic state peculiar to geometrical frustration. 

In this paper, we focus on the effects of the next- 
nearest neighbor (NNN) interaction on the pyrochlore- 
lattice classical Heisenberg model. Effects of the further 
neighbor interactions were investigated within a mean- 
field approximation by Reimers et al. up to the fourth 
neighbors [S]]- In the case of the NNN interaction J2 
only, they showed that a. q = order was stabilized for 
the antiferromagnetic NNN interaction J2 < 0, while an 
incommensurate q order was stabilized for the ferromag- 
netic NNN J2 > 0. Since then, several Monte Carlo (MC) 
simulations on this J1-J2 pyrochlore Heisenberg model 
were performed |9l-[ll|. 

For the antiferromagnetic J2 < 0, these works revealed 
that the system exhibited a first-order transition from the 
paramagnetic phase to the g = ordered phase with a 
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FIG. 1: (color online) A pyrochlore lattice. The nearest- 
neighbor interaction Ji and the next-nearest-neighbor inter- 
action J2 are indicated. 



collinear up- up- down- down spin structure [9|, |ll| . Selec- 
tion of such a collinear structure among possible q = 
states might be due to "order-by-disorder" mechanism, 
where thermal fluctuations lift the degeneracy through 
an entropic contribution to the free energy [2l|. 

For the ferromagnetic J2 > 0, on the other hand, 
the situation seemed to be more subtle. Tsuneishi et 
al. observed that the model with J2/J1 = —0.1 ex- 
hibited a first-order transition into a peculiar ordered 
state, where enhanced spin fluctuations apparently coex- 
isted with sharp Bragg peaks associated with a multiple- 
q structure Q • Chern et al. studied the case of smaller 
J2/IJ1I < 0.09 and showed that there occurred succes- 
sive phase transitions with an intermediate phase char- 
acterized by a finite nematic order parameter and a lay- 
ered structure, which was not predicted in the mean-field 



analysis [ll|. Chern et al. also pointed out that the low- 
temperature phase was a multiple-g state which were ba- 
sically the same as the one observed by Tsuneishi et al. In 
any case, the explicit spin configuration of the multiple-g 
state was not identified so far. 

In the present paper, we investigate the nature of the 
multiple-g state observed in the J1-J2 pyrochlore-lattice 
Heisenberg antiferromagnet with the ferromagnetic NNN 
interaction on the basis of a mean-field analysis and an 
extensive MC simulation. Particular attention is paid 
to the explicit spin configuration of the multiple-g state 
and the nature of spin fluctuations in such state. We see 
that there are mainly two stable multiple-g structures, 
the sextuple-g state and the quadruple-g state, each of 
which is characterized by whether they keep the cubic 
lattice symmetry or not. A first-order transition observed 
in earlier studies corresponds to the transition from the 
paramagnetic state to the cubic-symmetric multiple-g 
(sextuple-g) state. We predict that another phase transi- 
tion might occur at a lower temperature from the cubic- 
symmetric multiple-g state to the non-cubic multiple-g 
(quadruple-g) state. Refiecting their characteristic spin 
fiuctuations, internal fields in such multiple-g states ex- 
hibit a broad distribution. 

The rest of the paper is organized as follows. In Sec. 
im we describe our model and briefiy review the previ- 
ous results on the model. In Sec. IIIIl we present the 
result of our mean-field calculation determining the ex- 
plicit spin configurations of the possible multiple-q states. 
Our MC results are presented in Sec. IIVI The results are 
discussed in conjunction with the mean-field result, in- 
cluding the explicit spin configurations of the multiple-g 
states. Finally, we summarize our results in Sec. |Vl 



II. MODEL 

The model considered is the classical Heisenberg an- 
tiferromagnet on the pyrochlore lattice, whose Hamilto- 
nian is given by 



n = 



-Ji Yl ^^ 

(id) 



5,-J2 5] ^.-S,, (1) 



that the model exhibited a phase transition from the 
paramagnetic phase to a magnetically ordered phase, 
where the magnetic long range order was characterized 
by twelve incommensurate wavevectors of the FCC Bra- 
vais lattice: (g*,^*, 0),(g*, -g*, 0), {0,q* ,q*),{0,q*,-q*), 
{q*,0,q*),{-q*iO,q*) and their minus, with g* ~ || [g]. 
Since the pyrochlore lattice consists of four FCC sublat- 
tices, each wavevector has four independent eigenmodes 
and only a certain combination of them becomes unstable 
at the transition point. 

In fact, there are infinitely many ways to mix such 
twelve critical modes. Thus, just the determination of 
unstable modes is not enough to specify the explicit spin 
configuration of the ordered phase. Indeed, Reimers et 
at did not determine the ordered spin configuration, 
nor specified even whether it was a single-g state or a 
multiple-g state. 

Multiple-g state is generally incompatible with the 
fixed spin- length condition \Si\ — 1 so that it is not fa- 
vored in the classical spin system at low enough tem- 
peratures where the fixed spin-length condition needs to 
be observed. Meanwhile, the unstable critical mode of 
the present model entails the magnitudes of frozen spin 
moments differ from one sublattice to the other so that 
even a simple single-g state is incompatible with the fixed 
spin-length condition. It means that the multiple-g state 
might have a higher chance to be stabilized, particu- 
larly at moderate temperatures where thermal fluctua- 
tions play a role. 

As mentioned, a multiple-g ordered state character- 
ized by multiple Bragg peaks in the spin structure factor 
was reported in previous MC simulations on the model 
[3, |ll|. Interestingly, Tsuneishi et al. also reported 
that this multiple-g ordered state accompanies only small 
amount of spin freezing, and some sort of spin fluctua- 
tions apparently coexisted with sharp Bragg peaks. The 
detailed spin structure of the multiple-g state, however, 
has not been clarified so far. Especially, the origin of 
the spin fiuctuation reported by Tsuneishi et al. remains 
unsolved. 



where the first and the second sums are taken over all 
NN and NNN pairs Ji and J2 (see FigH]), respectively. 
We suppose the antiferromagnetic NN interaction Ji < 
and the ferromagnetic NNN interaction J2 > (J1-J2 
model), and a to be the length of the cubic unit cell. 
Note that the pyrochlore lattice, which can be viewed as 
an FCC lattice formed by regular tetrahedra of a fixed 
orientation, is a non-Bravais lattice. 

The ordering of this J1-J2 model was first studied by 
Reimers et al. based on a mean-field approximation. 
There, the full density matrix of the system was ap- 
proximated by a product of single spin density matri- 
ces with local effective fields which were determined so 
as to minimize the free energy. Reimers et al. found 



Chern et al. reported that yet another phase, a par- 
tially ordered collinear phase, might be realized between 
the paramagnetic phase and the multiple-q phase for suf- 
ficiently small values of IJ2/J1I < 0.09 ^. This par- 
tially ordered collinear phase is characterized by a finite 
nematic order parameter and a layered spin structure. 
Note that this phase is different from the multiple-q or- 
dered phase discussed above where the system retains a 
magnetic long range-order characterized by sharp Bragg 
peaks. Although the partially ordered collinear phase as 
discussed by Chern et al. is also interesting, we focus in 
the present paper on the nature of the multiple-g ordered 
state and of the associated fluctuations in the multiple-q 
ordered phase. 



III. MEAN-FIELD APPROXIMATION 

In studying the nature of the multiple-q ordered state, 
we begin by revisiting a mean-field analysis done ear- 
lier by Reimers et al. Q. Reimers et al. constructed a 
Landau-type free energy F of the pyrochlore Heisenberg 
antiferromagnet within a mean-field approximation up to 
the quartic order, 

F/Ns = -4rin47r 
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From the quadratic term of the free-energy expan- 
sion, one sees that the normal mode corresponding to 
the maximum eigenvalue A^. becomes unstable at Tc — 
■|Aq. where q* is the critical wavevector, leading to a 
phase transition to the ordered state characterized by 
the wavevector q* . When the maximum eigenvalue is 
degenerate as in the present case, the ordered-state spin 
configuration still remains largely undetermined at the 
quadratic level. In such a case, there are infinitely many 
ways of mixing $'. s with keeping the order parameter 
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where Bq is the order parameter corresponding to the 
Fourier magnetization of sublattice /i (/i = 1, 2, 3, 4) given 
by 

'^^''^ = ]^E'^i'^°^P(-9-n^"^ (3) 



constant. In order to specify the explicit ordered-state 
spin configuration, one needs to go to the quartic term. 
Reimers et al. made such an analysis only for the special 
case of a simple q — ordered state, whereas such an 
analysis has not been made for more general cases of an 
incommensurate ordered state, which is the target of our 
following analysis. 



where S^ is the spin at the site i belonging to sublattice 
/i, r,j is the position vector of that site and Ns is the 
number of spins belonging to sublattice /i. The sum X^ioi 
runs over all q^s which satisfy <7i -|- <72 + 93 + 94 = 0, and 
J^" is the Fourier transform of the exchange interaction 
between sublattices fi and v. The quadratic term can be 
diagonalized by a unitary matrix Uq with the eigenvalue 

Aq, 
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where i indicates each eigenmode. By transforming the 
order parameter to normal modes $' 



q Z^^q q- 
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Now, we wish to go beyond the analysis of Reimers et 
al. to derive the explicit ordered-state spin configuration 
of the J1-J2 model within a mean-field approximation. 
For the case of Ji < and J2 > of our interest, the 
maximum eigenvalue of Jq appears in a symmetric direc- 
tion q = {q, q, 0) and its cubic-symmetry counterparts. 
Along this direction, Jq has a form 



Jq=2 



with 
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the Landau free energy F is reduced to 
F/Ns = -4rin47r 
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j(i) = JlCOs| + 2J2, 

j'^-* — Ji COS — I- J2 I COS 1- COS -q 

4 V 4 4 

j(3) = Ji-f 2J2COS^. 



The eigenvalues of this matrix are calculated as 



(9) 



Ai 



-2Ji cos- — 4J2, 



A2 = -2Ji - 4J2COS-, 



A± = (Ji + 2J2)(cos| + l)±J(J: 
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with the corresponding eigenvectors given by 
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-A± + 2 ( Ji cos f + 2 J2) 
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The maximum engenvalue hes in the A_|_ branch at a 
wavevector q = q* ~ |^. Note that, since |a+| ^ 1 

generahy in the corresponding eigenvector [/+ , a single-g 
state associated with such an eigenvector cannot satisfy 
the fixed spin- length condition \Si\ = 1 for all sublattices. 
Hence, even a pure single-g state cannot be realized at 
T = as in multiple-g states, though it might be realized 
at finite temperatures due to thermal fluctuations. 

Due to the cubic symmetry of the pyrochlore lattice, 
A+(q*) is twelvefold degenerate. Hence, a variety of 
multiple-g states might be possible in principle. Even 
if one is to identify two eigenmodes corresponding to 
q and — q, there are still six independent eigenmodes. 
Within the Landau- type free-energy expansion ([B]), the 
free-energy difference among possible multiple-g states 
arises from the quartic term. 
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where the sum over i,j^k,l is taken over eigenmodes 
{i,j,k,l — 1,2,±). Just below the transition tempera- 
ture, one might neglect in the summation J2ijki S{q} ^he 
contributions from all modes other than the twelve criti- 
cal modes. Let us represent independent critical wavevec- 
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q^ = {0,q*,±q*) 

i±q*,0,q*)- (17) 

When one gives the number of critical modes to be mixed, 
there remain degrees of freedom associated with the am- 
plitudes and the directions of the critical modes *^±, 

which should be determined so as to minimize /j. 

In case of the single-q state, /j takes a minimum for a 
simple spiral state 

S{r) ex cos(q,f • r + e)ei ± sin(qrf • r + 6')e2, (18) 

where a = ± with ei _L 62 and |ei| = |e2| ~ 1, while the 
minimized f^ is given by 
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where a = a+(<7*). 

In addition to the simple single-g state, various types 
of multiple-g states are possible depending on the number 
of mixed critical modes. The multiple-g states with odd 
number of critical modes, i.e., the states with three or 
five mixed eigenmodes, turn out to have higher /4, and we 
consider here only the multiple-^ states with even number 
of critical modes, i.e., the states with two, four and six 
mixed eigenmodes, which are described as the double-g, 
the quadruple-g, and the sextuple-g states, respectively. 
For a given even number of mixed critical modes, one 
needs to minimize /j for a fixed order parameter m? with 
respect to the amplitudes and the directions of the critical 
modes $ + . Some of the details of this minimization 

procedure are given in appendix. 

In case of the double-9 state, the optimized spin con- 
figuration turns out to be a superposition of two dis- 
torted spirals of qf which have mutually orthogonal spi- 
ral planes as illustrated in Figl2] (a) (the dl state in ap- 
pendix). The minimized fi is calculated to be 
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(20) 
In case of the quadruple-g state, the optimized spin 
configuration turns out to be a superposition of two spi- 
rals q^ sharing a common spiral plane and two linearly- 
polarized spin-density waves (SDW), S{r) ex cos(q°' -r-t- 



(a) Double -q 



(a) Quadruple - 




(c) Sextuple -q 




FIG. 2: (color online) Schematic representation of several multiple-g states. Twelve arrows represent twelve wavevectors 
corresponding to critical eigenmodes of the model. A linearly-polarized spin-density wave (SDW) or a spiral depicted on each 
arrow demonstrates the type of the corresponding eigenmode. (a) The double-g state which is a superposition of two spirals of 
qf" with mutually orthogonal spiral planes, (b) The quadruple-g state which is a superposition of two spirals of qf" sharing a 
common spiral plane and two SDWs perpendicular to the spiral plane, (c) The sextuple-g state which is a superposition of six 
SDWs where each pair of q- shares a common axis, while axes of different qts are mutually orthogonal. 



0)e^, with its spin polarization perpendicular to the spi- 
ral plane as illustrated in Fig|2](b). The minimized fi is 
calculated to be 
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(21) 



Finally, in case of the sextuple-g state, the optimized 
spin configuration turns out to be a superposition of six 
linearly-polarized SDWs where each pair of q^ shares 
a common axis e^, and axes of different q^s are mutu- 
ally orthogonal ei _L 62 J- 63 as illustrated in Figl2] (c). 
Note that this sextuple-g state retains a cubic symmetry 
of the pyrochlore lattice, in sharp contrast to the other 
multiple-g states or the single-q state which break the 
cubic lattice symmetry. The minimized /4 is calculated 
to be 
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In FiglSl we show /4/r of various states discussed 
above, including the single-g state, the double-? state, 
the quadruple-g state and the sextuple-g state, as a func- 
tion of the parameter a. In fact, a depends on J2/J1 
only weakly. For the value of J2 < |Ji|, a takes values 
around a ~ 0.4: See the inset of Fig|3l Around this 
value, the free energy of the single-g state turns out to 
be much higher than those of the multiple-g states. This 
result is consistent with the observation of the multiple-g 
states in previous MC simulations [^, [ll|. Within the 
present mean-field approximation, the double-g state be- 
comes most stable for a ~ 0.4. The free energy difference 
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FIG. 3: (color online) Quadratic terms of the free energy /i 
of the single-5 state and of various types of multiple-g states, 
including the double-g state (the dl state in appendix), the 
quadruple-g state (the q2 state in appendix) and the sextuple- 
q state, are plotted versus the parameter a = a'^{q''): See 
Ea. p5|) . Inset exhibits a relation between a and J2/IJ1I. For 
a > 0.518 the double-g state and the quadruple-g state be- 
come locally unstable. Such unstable regions are indicated 
by thin dotted lines. For further details of the stability, see 
appendix. 



among various multiple-g states, however, turns out to be 
rather small, and fiuctuations neglected here might even- 
tually select the multiple-g state other than the double-g 
state. In order to determine the true ordered state of 
the model, we should carefully examine the effect of fluc- 
tuations. For this purpose, we perform extensive MC 
simulations of the model in the next section. 



IV. MONTE CARLO SIMULATION 

In this section, we present the results of our MC simu- 
lation on the Ji" J2 pyrochlore-lattice classical Heisenberg 
model with Ji < and J2 > 0. As a typical example, 
we deal with the case oi J2/J1 — —0.2. In the phase 
diagram reported by Chern et aZ.|llJ], this J2/\Ji\ ratio 
is sufficiently large so that the partially ordered coUinear 
phase should not appear. A direct transition from the 
paramagnetic phase to the multiple-g ordered phase is 
then expected. 

The pyrochlore lattice contains 16 spins in its cubic 
unit cell. The system size we deal with is of linear size 
L in units of cubic unit cell, i.e. the system contains 
N = 16L^ spins in total. Since the ordered state is gen- 
erally incommensurate with the underlying lattice, the 
system under periodic boundary conditions (BC) would 
be subject to severe finite size effects. In order to exam- 
ine such finite-size effects, we also consider free BC, and 
carefully compare the results between these two BC. In 
the case of free BC, we extend the lattice with a half of 
the unit cell in all three directions so that it keeps the 
cubic symmetry of the lattice. Then, the system contains 
TV = 2((2L + 1)3 + 1) spins in total. 

MC simulations are performed based on the standard 
heat-bath method combined with the over-relaxation 
method. The lattice size is of 8 < L < 24 for periodic 
BC, and of 8 < L < 32 for free BC. The system is grad- 
ually cooled from the high temperature. A run at each 
temperature contains typically (2 — 4) x 10^ MC steps per 
spin (MCS) and about a half of MCS are discarded for 
equilibration. Our 1 MCS consists of 1 heat-bath sweep 
and subsequent 10 over-relaxation sweeps. 

In FigHl we show the temperature and size dependence 
of the energy per spin. For both cases of periodic and free 
BC, an almost discontinuous change of the energy, a char- 
acteristics of a first-order transition, is observed. Such a 
first-order nature of the transition is consistent with the 
results of previous calculations [^, [ll| ■ We estimate the 
bulk transition temperature to be T^ ~ 0.178| Ji|. 

At temperatures below Tc, several types of metastable 
states are observed in our simulations, depending on 
the spin initial conditions and the random-number se- 
quences. In order to probe the spin configurations in 
these metastable states, we calculate the spin structure 
factor F{q) defined by 
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where (• • ■ ) represents a thermal average. From the spin 
structure factor, we find that none of metastable states is 
a single-g state. All of them are multiple-g states where 
the multiple magnetic Bragg peaks coexist in the spin 
structure factor. 

In FigjS] we show the typical spin structure factors 
of these metastable states calculated at a temperature 
T = 0.17|Ji| <Tc~ 0.1781 Jil under periodic BC. These 




FIG. 4; (color online) Energy per spin versus the tempera- 
ture for both cases of periodic boundary conditions (a), and 
free boundary conditions (b). The interaction parameter is 
J2/J1 = -0.2. 



metastable states might be classified into two types ac- 
cording as whether they keep the cubic symmetry or 
not. In the cubic-symmetric state, there are six inde- 
pendent main Bragg peaks with the same intensity at 
q = —{h*, h* , 0) and at its cubic-symmetry counterparts 
with h* ~ 5/4 (see FigjS] (a)). The peak position h* is 
close to the corresponding mean- field value h* ~ 1.263. 
There are also other Bragg peaks related to the main 
peaks via the reciprocal lattice vectors of the FCC lat- 
tice. 

Since the observed cubic-symmetric state involves six 
critical wavevectors, we identify this state as a sextuple- 
q state. In case of periodic BC, we sometimes observe 
"almost cubic symmetric" states where the position of 
the Bragg peak shifts by ^ , or one of the Bragg peaks 
splits into two, due to the finite size effects associated 
with periodic BC. We regard them as modified forms of 
the sextuple-g state. 

In the observed non-cubic state, two out of six main 
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FIG. 5: (color online) The spin structure factor of the multiple-g states along the directions of q = —{h,h,0) and of its 
cubic-symmetry counterparts, calculated at T = 0.17|Ji| with J2/J1 — —0.2. The lattice is of L = 16 with periodic boundary 
conditions. Each curve corresponds to different directions of wavevectors related via the cubic symmetry of the lattice, (a) 
The sextuple-g state where Bragg peaks appear in all six directions with the same intensity, (b) The quadruple-q state (of 
type 1) where Brag peaks appear only four out of six directions. Half of main Bragg peaks corresponding to the wavevectors 
in the same plane (/i, ±/i, 0) are larger in magnitude than those in the other plane (0, ft, ±/i). (c) The quadruple-g state (of 
type 2) where Brag peaks appear in four out of six directions with the same intensity. In this type 2 structure, all Bragg peaks 
including main and submain are split into two, one at {q,q,0) and the other at (q — ^,g, 0). The peaks shown in the figure 
are the one at (g, q, 0). 
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FIG. 6: (color online) The spin structure factor of the multiple-g states along the directions of q = —{h,h,0) and of its 
cubic-symmetry counterparts, calculated at T = 0.17|Ji| with J2/J1 = —0.2. The lattice is of L = 32 with free boundary 
conditions. Each curve corresponds to different directions of wavevectors related via the cubic symmetry of the lattice, (a) 
The sextuple-g state where Bragg peaks appear in all six directions with the same intensity, (b) The quadruple-g state (of 
type 1) where Brag peaks appear only four out of six directions. In (b), the main and sub peaks appear along the direction 
q = —{h,h + 6,0), slightly off the high symmetric direction {h,h,0), where S — 1/{L + 1) is simply the mesh size of our 
measurements in the wavevector space. The peak height shown here is of the one at (h, h + S,0), while the peak position is set 
at /i = ^/[h+{h + Sy]/2. 



Bragg peaks with wavevectors on a common plane, e.g., 
{q,0,q) and {—q,0,q), vanish, while other four remain. 
Since there are four critical wavevectors in this non-cubic 
state, we identify the state as a quadruple-g state. Two 
different types of the non-cubic states are observed, de- 
pending on the intensity ratio among the remaining four 
main Bragg peaks. In one type (type 1), two out of four 
remaining Bragg peaks with wavevectors on a common 
plane, e.g., {q,q,0) and {q,—q,0), have stronger inten- 
sity than the other two with wavevectors lying in the 
other plane, e.g., {0,q,q) and {0,q,—q): See Figl5] (b). 
In the other type (type 2), four main Bragg peaks have 



equal intensities (see FiglS] (c)). In this type 2 struc- 
ture, all Bragg peaks including main and submain are 
actually split into two, one at (g, q, 0) and the other at 
{q-^,q,0). 

In Figini we show the typical spin structure factors 
calculated at the same temperature T = 0.17|Ji| under 
free BC. Similarly to the case of periodic BC, we find a 
cubic-symmetric state (Fig [5] (a)) and a non-cubic state 
of type 1 (Fig|6] (b)). By contrast, the non-cubic state 
of type 2 is never observed in contrast to the case of 
periodic BC. (In the non-cubic state observed under free 
BC, the main and sub peaks appear along the direction 
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FIG. 7: (color online) The Monte-Carlo time evolution of 
the cubic-symmetry-breaking order parameter nic defined by 
Eq.dMll, calculated at T = 0.15|Ji| for a L = 48 lattice 
with free boundary conditions. The interaction parameter is 
J2/J1 = —0.2. Initially, a half side of the system is prepared 
to be a sextuple-g state, while the other half a quadruple-g 
state. The two curves represent the subsequent time evo- 
lution of the cubic-symmetry-breaking order parameter rric, 
each calculated at each half side of the system. The data in- 
dicates that the cubic state expands as the simulation goes 
on. 
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FIG. 8: (color online) Energy difference between the sextuple- 
q state and the quadruple-g state plotted versus the tempera- 
ture for both cases of periodic and free boundary conditions. 
The interaction parameter is J2/J1 ~ —0.2. In case of pe- 
riodic boundary conditions, we choose a pair of metastable 
states which have the lowest energy among several "almost 
cubic-symmetric" states or "almost quadruple-g (type 1)" 
states where peak shift or peak splitting occur due to finite 
size effects. 



of q = —{h,h + S,0) slightly off the high symmetric 
direction {h, h, 0) reflecting the non-cubic character of 
the ordered states. However, the magnitude of the shift 
d is always of order of 1/L, apparently vanishing in the 
thermodynamic limit.) 

We deduce that the type-1 state rather than the type- 
2 state is a stable quadruple-g state due to the follow- 
ing two reasons. First, the type-1 quadruple-g state ap- 
pear as a metastable state in both periodic and free BC, 
while the type-2 quadruple-g state is never realized in free 
BC. Second, within our mean-field analysis, the type-1- 
like state with the spin structure factor of unequal peak 
heights, which corresponds to the q2 state in appendix, 
has a lower free energy than the type-2-like state with 
the spin structure factor of equal peak heights, which 
corresponds to the ql state in appendix. 

In contrast to the mean-field result for J2/J1 — —0.2, 
we do not observe in our MC simulation the double- 
g metastable state below Tc- Since the double-g state 
can continuously be changed into the quadruple-g state, 
the double-g state is probably unstable toward the 
quadruple-g state in this parameter region. 

Next, we wish to determine either the cubic sextuple-g 
state or the non-cubic quadruple-g state is stable below 
Tc- Since both states are metastable states surviving for 
a very long time in the course of MC simulation, it is 
extremely difficult to identify which state is thermody- 
namically stable by performing fully thcrmalized simu- 
lation below Tc- In order to identify the truly stable 
state below Tc, we then employ the mixed-phase method 
[111 |23 | . In this method, one prepares as an initial state 
a two-phase coexisting state, where the sextuple-g state 
or the quadruple-g state occupies each half of the lattice. 
By thermalizing such a state and by monitoring which or- 
dered state expands during the course of the subsequent 
simulation, one can determine which phase is more sta- 
ble. 

In Figl?! we show an example of the MC time evolution 
of such a two-phase coexisting initial state in our MC 
simulations. We look at the time evolution of the cubic- 
symmetry-breaking order parameter rric for each half of 
the lattice, defined by 

^ ^ (max {Y,i S., ■ Sj+s)} - (min {J2, S., ■ Sj+s)} 

" ~ (max {J2i St ■ Si+s)) + (min (^^ S^ ■ Si+s)) ' 

(24) 
where max {J2t Si ■ St+s) and min {J2i Si ■ Si+s) repre- 
sent the maximum and the minimum values of the sum of 
the inner product Si ■ Si+s among three nearest-neighbor 
bonds, denoted by S, emanating from a given site i. 
This quantity measures the extent of the cubic-symmetry 
breaking in the spin configuration: It tends to be large 
for a non-cubic quadruple-g state, while it tends to be 
small for a cubic sextuple-g state. As can be seen from 
the MC time dependence of rric shown from FigEl at 
the measuring temperature T — 0.15|Ji|, the sextuple- 
g state expands its area, and the system finally settles 
in the sextuple-g state with smaller ijic. By performing 
such a mixed-phase method for various system size up to 



L = 24 for periodic BC, and up to L = 48 for free BC, we 
find that the sextuple-q state is stable at least in the tem- 
perature range T > 0.15Ji. Unfortunately, the dynamics 
of the system becomes so slow with further lowering the 
temperature that the domain wall between the two or- 
dered states hardly moves, which hampers to determine 
the stable state at lower temperatures T < 0.15Ji. 

In the low temperature limit, the entropy contribution 
to the free energy becomes smaller, while the energy con- 
tribution becomes dominant. Thus, the state which has 
the lower energy should be realized. In FiglH we show the 
energy difference between the sextuple-g state and the 
quadruple-g state versus the temperature. For both cases 
of periodic and free BCs, the energy of the quadruple-g 
state becomes lower than that of the sextuple-q state at 
low enough temperatures. It means that the quadruple-g 
state should be realized as a stable state at sufficiently 
low temperatures. 

By combining the results of the energy analysis with 
those of the mixed-phase method, we conclude the or- 
dering behavior of the J1-J2 model with J2/J1 = —0.2 
as follows. With decreasing the temperature, the system 
first exhibits a phase transition from the paramagnetic 
phase to the sextuple-g ordered state at a temperature 
T = Tci ~ 0.178|Ji|. The transition is of first order. 
Note that the cubic symmetry of the lattice is still fully 
respected even below Td in spite of the existence of the 
magnetic long-range order characterized by sharp Bragg 
peaks. With further decreasing the temperature, another 
phase transition occurs aX T — Tc2 < 0.15|Ji| into the 
quadruple-g state. It accompanies the breaking of the 
cubic symmetry of the lattice. This transition is proba- 
bly of weakly first order. 

Now, we turn to the determination of the detailed spin 
configurations in the multiple-g states identified above. 
We begin with the cubic-symmetric sextuple-g state. 
Guided by our mean-field result and by inspecting the 
spin configurations obtained by MC simulations, we ob- 
serve that the ordered spin configuration in the sextuple- 
q state can be represented by the superposition of six 
linearly-polarized SDWs as 

Si^^r) =l\u\ cos{q+ • r + 61+) 



U^_ cos{q^ ■ r + 9^) 



+ C/^_ cos(q7 • f 
-I- [/'"_ cos(q3" • r 

^3 



(25) 



where S'-'^^r) = (si^^r) , S^^^r) , si'"\r)) represent the 
spin at the position r belonging to sublattice fi. Any 
spin configuration generated from the one given by (j25p 
via global spin rotation in spin space is equally allowed. 



The coefficient Uj^ takes a value unity or —a, depending 
on whether the site on the sublattice /i has a NN bond 
in the direction of the wavevector q ~ qf [U^^ = 1) or 

not (C/'^± = —a). Indeed, the spin configuration given 

by Eg. ([25]) is fully consistent with the mean-field result 
shown in Sec. Ill and appendix. 

In the the case of the quadruple-g state, inspection of 
the MC data lead us to the ordered-state spin configura- 
tion represented by the superposition of two spirals and 
two linearly-polarized SDWs as 



Si'Hr)^Ixy 



C/^cos(q+•r + 0+) 
+ U^^cosiq^-r + d^) 

U^. sin(q+ -r + et) 
+ C/^^_sin(gr.r + 0r) 
il^. cos{q+ ■ r + e+) 

"2 



(26) 



where the two wavevectors with stronger intensity in the 
spin structure factor are associated with q^ and those 
with weaker intensity with q^- Any spin configuration 
generated from the one given by (j26p via global spin ro- 
tation in spin space is equally allowed. The same rule 
given above for the sextuple-g state is understood also 
for the coefficient U^. This spin configuration is also 
fully consistent with the mean-field result of Sec. Ill and 
appendix. 

In order to investigate the nature of fluctuation in the 
ordered state, we compute the time-dependent spin over- 
lap g(^^ [t) defined by 



g(2)(i)^ 



Q,/3 



EE^i^^(-.. 



(m) 



io)4^V., 



M 



to + t) 



(27) 

which measures an overlap between the spin configura- 
tions at time to and at time to -I- t. Note that q^'^\t) is 
defined so as to be invariant under any global spin rota- 
tion in spin space. 

We also employ this quantity to check the validity of 
the proposed spin structures given by Eas. ([25l) and (|26p . 
For this purpose, we compute g"' (t) in the following two 
ways: The two spin configurations at time to and at sub- 
sequent time to + t are taken, either (i) from the raw spin 
configuration data of our equilibrium MC simulations, 
or (ii) from the the proposed model spin structure de- 
scribed by Eq. (P51) or (^5)1 which is evaluated by Fourier 
transforming the raw MC data and by extracting the 
amplitudes (/, I^y and Iz) and the phases (0^) associ- 
ated with the critical mode at g^ . In the procedure (ii), 
the a-factor associated with U or U is also taken as a 
fitting parameter, not a given constant, while the contri- 
bution of other non-critical modes are simply neglected. 
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FIG. 9: (color online) The Monte-Carlo time dependence of the overlap q (t) defined by Eq. (pTJ) . calculated at T = 0.17|Ji| 
for a 1/ = 16 lattice with periodic boundary conditions, in the cases of (a) the sextuple-g state and of (b) the quadruple-q state. 
The interaction parameter is J2/J1 ~ —0.2. Two curves in the figure represent q (t) calculated either directly from the raw 
Monte Carlo data or from the proposed model spin structure evaluated from the Monte Carlo data: See the text for further 
details. In both cases of (a) and (b), the two curves agree with high precision. 



If the proposed model structures properly represent the 
ordered-state spin configuration, q^^\t)s calculated in the 
above two ways (i) and (ii) should agree. 

In Fig[Sl we show the MC-time dependence of q^'^\t) 
calculated in the two ways (i) and (ii) for the sextuple-q 
state (a), and for the quadruple-g state (b). Indeed, in 
both cases of the sextuple-g and the quadruple-g states, 
g(^) (i) calculated in the two ways agree with high preci- 
sion, indicating that our model spin structures (|25p and 
((26)) properly represent the actual ordered-state spin con- 
figurations. Second, the computed q^'^^t) exhibit a sig- 
nificant time variation indicating that the spin config- 
uration changes considerably in the course of our MC 
simulation. Indeed, as can be seen from FigfTOl this time 
dependence comes from the time dependence of the phase 
factor 9- , whereas the amplitude turns out to be nearly 
time- independent . 

Establishing the spin configuration of the ordered 
state, we wish to further investigate the nature of spin 
fluctuations in the ordered state. Two kinds of spin fluc- 
tuations associated with the identified spin structures are 
observed. The first one is a fluctuation of the phase fac- 
tors 9f^ as mentioned above. While the amplitudes /, 
I'xy , Iz , which contribute to the Bragg intensity, turn out 
to be nearly constant in time (see FigHU] (a) and (c)), 
the phase factors Of^ fluctuate a lot even when the sys- 
tem exhibits sharp Bragg peaks (see Fig[TU](b) and (d)). 
The second type of spin fluctuation is a necessary out- 
come of the distribution of ordered spin moments in the 
multiple-g state, where some spins should reduce their 
frozen moments exhibiting large fluctuations. 

First, we examine the effect of phase fluctuations. As 
we have seen in the dynamics of 0^{t) shown in FigfTUl 
(b) and (d) , phase fluctuations do occur in finite-size sys- 



tems. The question to be addressed is whether such 
phase fluctuations disappear or not in the thermody- 
namic limit. For this purpose, we define the phase auto- 
correlation function by 



Cit) = {cosidit)~9{0))), 



(28) 



where 9{t) represents Of{t). In Fig[TT](a), we show the 
phase autocorrelation function of 61 in the quadruple- 
q state. As can be seen from the figure, the relaxation 
of the phase becomes slower with increasing the system 
size, suggesting that the phase motion might eventually 
be locked in the thermodynamic limit. The computed 
phase relaxation follows a simple exponential decay char- 
acterized by the phase correlation time t. In FigfTTKb). 
we plot 1/r as a function of the inverse system size 1/N . 
As can be seen from the figure, the phase relaxation time 
diverges almost linearly with iV, and we conclude that 
the phase motion will eventually be locked in the ther- 
modynamic limit. Hence, spin fluctuations arising from 
phase fluctuations are finite-size effects and are expected 
to vanish in the thermodynamic limit. 

Next, we discuss the second source of spin fluctua- 
tions. Even when phase fluctuations vanish in the ther- 
modynamic limit, spin fluctuations associated with the 
multiple-q nature of the ordering should still remain. Re- 
member that the multiple-g order entails the the exis- 
tence of certain fraction of spins with reduced frozen mo- 
ments, which should be caused by significant local spin 
fluctuations. 

Such spin fluctuations intrinsic to the multiple-q or- 
der might be observable experimentally via local probes 
like NMR. In FigfT^ we exhibit the distribution of inter- 
nal flelds in the sextuple-g (Eas.(|^5|)) and the quadruple- 
q (Eq. (|^5|) ) states probed at an O site of a pyrochlore 
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FIG. 10: (color online) The Monte-Carlo time dependence of the amplitudes I,Ixy,Iz ((a) and (c)) and the phase factors 6- 
((b) and (d)) for the cases of the sextuple-g state ((a) and (b)) and of the quadruple-g state ((c) and (d)), as obtained by 
fitting the raw Monte-Carlo data to the model spin structures of Eas. (|25|l and p6|) . Concerning the amplitudes, plotted here 

are the average intensities over sublattices fi, i.e., /av = \/X]u [^'^q] /4- The temperature is T = 0.17|Ji| and the lattice is 
L = 16 with periodic boundary conditions. The interaction parameter is J2/J1 = —0.2. The amplitudes, which contribute to 
the Bragg intensity, are kept nearly constant in time, whereas the phase factors fluctuate a lot. 



oxide A2B2O6O' with the oxygen-location x parameter 
X = 0.319 [23]. We assume here that the B site is mag- 
netic, and internal fields are borne by the dipolar interac- 
tion. As can be seen from the figure, both the sextuple-q 
and the quadruple-q states exhibit broad internal-field 
distributions, in contrast to the one associated with a 
simple all-in/all-out structure shown in the inset. The 
latter exhibits a sharp delta-function-like distribution. 
The observed broad distributions of internal fields mean 
that the magnitude of frozen spin moments are spatially 
distributed, as expected from the multiple-g nature of the 
ordering. 



V. SUMMARY AND DISCUSSION 

We studied the nature of the multiple-^ ordered states 
realized in the J1-J2 pyrochlore-lattice Heisenberg model 
by means of a mean-field analysis and a Monte Carlo 
simulation. 

Performing a mean-field analysis beyond the previous 
analysis by Reimers et aL[^, we could explicitly deter- 
mine the possible multiple-q spin structures of the model. 



which include a cubic-symmetric sextuple-g state, a non- 
cubic quadruple-g state and a non-cubic double-g state. 
We have found that these multiple-g states have consid- 
erably lower free energy than that of the single-g spiral, 
while the free energy difference between these different 
multiple-g states are rather small. 

With reference to the mean-field results, we also per- 
formed an extensive MC simulations of the model, mainly 
for the case of J2/J1 = —0.2, to determine which state is 
really stabilized as an ordered phase. As expected from 
the mean-field analysis, we found that the system exhib- 
ited a phase transition from the paramagnetic phase to 
the multiple-^ ordered state characterized by the multiple 
peaks in the spin structure factor. Recent studies have 
revealed that such multiple-g states are also stabilized in 
other frustrated Heisenberg magnets, e.g. the triangular 
lattice Heisenberg antiferromagnet with the next-nearest 
or the third neighbor interactions under magnetic fields 
[24| . With the help of the mixed-phase method, we found 
that the cubic-symmetric sextuple-g state is stabilized 
just blow the transition temperature T = Td ~ 0.178| Ji| 
down to at least T = 0.15|Ji|, while, at sufficiently 
low temperatures, the non-cubic quadruple-g state be- 
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FIG. 11: (color online) (a) The Monte-Carlo time dependence 
of the phase df" autocorrelation function in the quadruple-g 
state at T = 0.17Ji . The lattices are L = 8, 12, 16 and 20 with 
periodic boundary conditions, (b) The phase autocorrelation 
time obtained as in (a) is plotted versus the inverse system 
size 1/A'^. Dotted line exhibits a liner fit of the data. 



comes stable. Hence, another phase transition from the 
sextuple-g state to the quadruple-g state should occur at 
a temperature T — Tc2 < 0.15| Ji|. 

The nature of each transition was also examined. The 
transition at T = Td is first order, as was already indi- 
cated by previous studies. Note that, at this transition 
the cubic symmetry of the lattice is still fully preserved. 
It is remarkable that, in spite of the incommensurate and 
rather complex nature of the spin order, the cubic sym- 
metry of the lattice is fully respected in the ordered state. 
The second transition at T = Tc2 is weakly first order, 
and it accompanies a spontaneous breaking of the cubic 
symmetry of the lattice. 

It should be noticed that the phase transition between 
the sextuple-g and the quadruple-q states analyzed in the 
present paper is distinct from the transition between the 
paramagnetic state and the nematic state (or the one be- 
tween the nematic state and the multiple-g state) as dis- 
cussed by Chern et a?.[ll|. We also performed MC sim- 




FIG. 12: (color online) An internal field distribution at an O 
site of a pyrochlore oxide A2B2O6O' with the oxygen-location 
X parameter x = 0.319 where B site is magnetic, calculated 
for an ideal sextuple-q state (|25[) and quadruple-g state (|26|) . 
respectively. Dipolar interactions are assumed as an origin 
of internal fields. Inset exhibits an internal field distribution 
corresponding to a commensurate all-in/all-out structure on 
the pyrochlore lattice. 



ulations for several smaller values of J2/I Ji| < 0.1 where 
Chern et al. observed the nematic phase. In such smaller- 
J2/\Ji\ region, we found as metastable states both the 
sextuple-g and the quadruple-g states in the parameter 
region where Chern et al. reported the multiple-^/ state 
(Chern et al. did not specify the type of the multiple-g 
state). Thus, there remains a possibility that, even for 
smaller J2/IJ1I- values where the nematic phase appears 
at higher temperature region, another phase transition 
from the sextuple-g state to the quadruplc-g state occurs 
at a lower temperature. 

We also determined the explicit spin configuration of 
the multiple-g ordered state. The sextuple-q state is a 
superposition of six SDWs running along the wavevector 
{q*,q*,0) and its cubic-symmetry counterparts. In this 
sextuple-g state, the system fully retains a cubic sym- 
metry of the lattice. By contrast, the quadruple-g state 
is a superposition of two helices and two SDWs. This 
state spontaneously breaks the cubic symmetry of the 
lattice as one can easily see from the fact that two out of 
six critical wavevectors should vanish in this state. Both 
the sextuple-(7 and the quadruple-g spin configurations 
are consistent with the ones obtained in our mean-field 
analysis. 

In both cases of the sextuple-g and the quadruple-^ 
states, a broad distribution of internal fields was ob- 
served. Such broad distributions reflect the distribution 
of frozen spin moments which arises from the multiple-g 
nature of the ordering. Hence, if the pyrochlore mag- 
nets are in the multiple-q ordered phase we proposed, one 
should observe two apparently conflicting features, i.e., 
the coexistence of sharp Bragg peaks measured by neu- 
tron diffraction and of enhanced spin fluctuations mea- 
sured by local probes such as NMR. Experimental obser- 
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vation of the multiple-g states as revealed here remains 
most interestmg. 
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Appendix: details of the mean-field approximation 

In this appendix, we present the details of our mean- 
field analysis. We consider here twelve critical eigenvec- 
tors only, corresponding to six wavevectors of ([T7|. As 
mentioned in section III, the free-energy difference be- 
tween different types of multiple-g states arises at the 
quartic term of the free energy P^. Thus, we calculate 
here the quartic term of the free energy f^ for several 
multiple-g structures. 

For simplicity, we abbreviate ^'I^ as ^ia-, where i — 
1,2,3 and a — zt. The quartic term /4 consists of follow- 
ing four terms, Ai, A2, A'2 as and A4 as 



/4- 



9T 



80(l-(-a2)' 



(2 + 2q*)Ai -f 32q2A2 



^{l + 2a^ + a'^)A'2 + Ma^Ai , (A.l) 



where Ai may be regarded to represent the interaction 
among i distinct wavevectors, and are given by 
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4$4 



(A.2) 
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+ 1*^^ • *j„ 
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. ■ **.' ' , 



(i^j) a,a' 



(A.4) 



and 



Ai^ Y^ 'Re 






(A.5) 



where ^ia = |*icr| and Re[---] means a real part, 
while J2(i^j) means the summation over {i,j) = 
(1,3), (2, 1), (3, 2). The A4 term represents interactions 
of four wavevectors such as q'^ , qf , q^ ,q^ , which satisfy 



Qt + Qi -it +Q3 = 0- 



1. The single-g state 

In the case of the single-g state, only the Ai term con- 
tributes to fi. If we fix the amplitude <I>?g. = m^ associ- 
ated with the wavevector qf , f^ takes a minimum when 
^ia satisfy 

*,, • *,, = 0. (A.6) 

This condition means that ^ia- is given by 

*,„ = me^''(ei+ie2), (A.7) 

where ei and 62 are orthogonal unit vectors in spin space. 
In the real space, this ^ia means a spiral, 

S{r) (X cos(qf • r + e)ei - sm{qf ■ r + 9)e2. (A.8) 



By substituting this condition into Eq. ljA.ip . we get 
the minimized free energy for the single-q state as 



p(si) 



gr 



4\™4 



fl'''> = ^(l + a^)m 

' 10(l + a2)'^ ^ 



(A.9) 



2. The double-q state 

In mixing two wavevectors, there are two possible 
ways: One is to mix a pair of q^ and q~ , and the other 
is to mix <7j and q- with i 7^ j. For a < 1, the former 
always minimizes the free energy, since the coefficients of 
A2 is smaller than that of A'2 . Hence, we consider below 
a mixture of q^ and q^ . In this situation, the yl2-term 
vanishes identically. 

The j42-term is minimized for 



*»+ *j 



0, 
*.+ • **_ - 0. 



(A.IO) 



If we assume that $i+ forms a spiral given by Ea. (IA.7l) . 
these conditions are satisfied when $i_ is perpendicular 
to the spiral plane formed by $i+, i.e., 



*i_ ex e^'ea. 



(A.ll) 



where 63 is a unit vector perpendicular to the spiral plane 
formed by #i+,ie., e^ — ei x 62. In this situation, the 
spin configuration corresponding to $j_ forms a spin- 
density wave (SDW) 



S'(r) ex cos(qj ■r + 9)e3. 



(A.12) 



Hence, the A2-term favors a mixture of a spiral and a 
SDW. 
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As we have seen for the case of the single-g state, Ai is 
minimized when ^i^- forms a spiral. The competition be- 
tween the Ai- and the A2-terms often causes a distortion 
of a spiral structure, i.e., 






(Ri, 



'^-'■icr ) : 



where Ria and lia are orthogonal real vectors 



■^irr -L ^h 



(A.13) 



(A.14) 



while the amplitudes Ria — \Ria\ and Ua — \Iia\ are 
generally different from each other. Note that such a 
distorted (or elliptical) spiral can continuously be trans- 
formed to an isotropic spiral {Ria = ha) and a SDW 
[Ricr = or lia = 0). Thus, we assume both $i+ and 
#i_ form distorted spirals Eq. (jA.13p . and minimize /4 
with respect to the directions of spiral planes, the phase 
factors 9ia, and the amplitudes Ria and lia. 



Under this situation, the A2-term is minimized when 
the two spiral planes are mutually orthogonal, i.e., 



(i?,+ X /,+ ) ± (i?,_ X /,_). 



Then, $i+ and $i_ are given by 



(A.15) 



*,+ = e*''"+ [Ri+ei + iKi+Ii+e2) , 

*,_ = e'"^- (i?,_e3 + zK,_/,_ei) , (A.16) 

where Kia = ±1 represents the chirality of each spiral. 
In this structure, the two spirals share one of the basis 
vectors ei while the other two basis vectors 62 and 63 
are mutually orthogonal. 

By substituting ^ia into Eq. ()A.l[) . /4 is obtained as 



9T 



80(l + a2)-^ 



2(1 



Mm 



$: 



i+ 



$? 



i?^ - Ilf + {Rl Ilf] } + 32a' \<i>l^l + 2i?^it] 



(A.17) 



Note that this /4 is independent of the phase factor 9ia 
and the chirality Kia. By minimizing f^ with keeping 
$f_i_ -I- <i>|_ — m? fixed, we find that there are two types 
of double-g structure, i.e., a mixture of two spirals (the 
dl state) and a mixture of a spiral and a SDW (the d2 
state). 

In the case of a mixture of two spirals (the dl state), 
/4 takes an extremum 



(•(dl) - 



9r 



for 



40(H-a2)' 



Ri+ 

h+ 

Ri-i- 



3 + 4a^ 



3a^- 



1 



(l + a^) 



(A.18) 



= Ri- , 



1 + a^ 



4(l + a2) 

1 + 4a2 -I- 

4(l + a2)2 



^m'. 



(A.19) 



In the dl state, two spirals are distorted {Ri+ < Ii+, 
h- < Rl-) to lower the A2-term. For a^ < a^ == (2-\/3) 
(ttc — 0.518), this dl state becomes a minimum of f^, 
while for a > etc it becomes locally unstable. 

In the case of a mixture of a spiral with q^ and a SDW 
with q~ (the d2 state), f^ takes an extremum when the 
spiral is isotropic and the SDW is orthogonal to the spiral 



as 



Ri+ — Ii+, 
h- = 0, 



R 



3 - 4^2 -h 3a^ 

'i 

2 (5 - 8a2 + 5a4) 
,2 _2-Aa' + 2a'^ ^ 



^'- 5 - 8a2 + 5q,4 



(A.20) 



where fi is calculated to be 



„(d2) _ 
/4 — 



9T 



20(1+^2)2 



(3 - 4a2 + 3^4)2 
5 - 8a2 + 50-4 

+ 3(l + a4) 



(A.21) 



For a2 > ^2 = (2 - V3) {a^ ~ 0.518), this d2 state 
becomes a minimum of f^, while for a < ac it becomes 
locally unstable. 

In Fig[T31 we plot /4/T of the single-g spiral state (the 
s state) and of the two kinds of double-g states (the dl 
and d2 states) as a function of a. Note that the critical 
value ac {(x^ = 2 — v3) of the stability is in common 
between the dl state and the d2 state. Thus, the mixture 
of two spirals (the dl state) is stable for a < ac — 0.518, 
while the mixture of a spiral and a SDW (the d2 state) 
is stable for a > ac. Since a is in the range of a ~ 0.4 
for the realistic parameter value of J2/ Ji, the mean- field 
approximation favors the dl state over the d2 state in 
the J1-J2 model. 



15 



0.58 




\/i^' 


- 


0.56 




\ 


- 


^0.54 

£:!. 0.52 


" " - . 


\ 




< 0.5 


/■(q2) y^ 

J A yy 

yy Ml) 


\S 


- 


0.48 


j(qi) 

/4 - 


0.46 


,(d2)\ 
J4 


V. 



0.2 



0.4 



0.6 



a 



FIG. 13: (color online) Quartic terms of the free energy /4 
of the single-g spiral state (the s state), of the two kinds of 
double-g states (the dl and d2 states) and of the two kinds of 
quadruple-g states (the ql and q2 states) are plotted versus 
the parameter a. For a > 0.518, the dl state and q2 state 
become unstable, while for a < 0.518 the d2 and the ql states 
become unstable. The unstable regions are indicated by the 
thin dotted lines. 



3. The quadruple-g state 

When four wavevectors coexist, the A4-terin con- 
tributes to the free energy. Since the increase of the 
number of mixed wavevectors tends to enhance the con- 
tribution of the A2- and the j42-terms, it is necessary to 
lower the ^4-term to stabilize the quadruple-g state. As 
an example, we consider here a mixture of Qi and q^ . 

Since spiral structures tend to lower the Ai-term while 
SDW structures tend to lower the A2- and the Aj-terms, 
we again assume that four $io-s form distorted spirals 
described by Eq. (jA.13[) , and minimize f^ with respect to 
the directions of spiral planes, the phase factors ^io-, and 
the amplitudes Ria- and Iia- 

First, we optimize the spiral plane. In order to lower 
the A2- and the A2-terms, the ideal situation would be 
that four spiral planes are mutually orthogonal. However, 
this is simply impossible for the Heisenberg spin. Since 
the coefficient of A2 is smaller than that of A2, next way 
might be that a pair of qf (or q^) shares a common 
spiral axis which is orthogonal to the others like. 



{Ri+ X Ji+) II (ili_ X /i_), 

{Ri+ X /i+) ± (^3+ X J3+), 
iRi+ X Ii+) ± {R3- x/3_). 



(A.22) 



Remaining spiral planes of qf prefer to be orthogonal to 
each other if we consider only the yl2-term. However, if 
we consider also the A^-tcTin, the sum of the A2- and 
the A4-terms takes a minimum when the pair of q^ also 
shares a common spiral plane orthogonal to that of q^ 
Uke 



(i?3+ X /3+) II {R3- X I3 



(A.23) 



as long as they are not reduced to the double-g or the 
single-q states. Thus, four ^i^s are given by 



*3 



= 6'"'+ {Ri+ei + iKi+h+e2) , 
= e}^^- [Ri-ei + iKi_/i_e2) , 



o«f3 + 



(Ra-Bs + iK3_/3_ei^ 



(A.24) 



where Kia- = ±1 represents the chirality of the spiral. 

Next, we optimize the phase factor, 9i^, and the chi- 
rality, KiCT- Although the Ai-, the A2-, and the A'2-tcTms 
are independent of the phase factor and the chirality, the 
v44-term depends on them. By substituting Ea. (|A.24p 
into Ea. (|A.5p . the y44-term is calculated as 



At = cos (6*1+ + 6li_ - 6I3+ + 6l3_) 

(i?l + i?l_ - Kl + Kl_/l+/l_) 
X {R3+R3- + K3+K3_/3+/3_) 

4- 2K3+K3_ii:i + i?i_/3+/3_ 

Then, the ^4-term is minimized for 



(A.25) 



cos (t^n 



where A4 is given by 



9.3+ +03-) = -!, 
Kl+Kl- = -1, 
K3+K3- = 1, 



(A.26) 



A4 = — (i?i+i?i_ + /i+/i_) (i?3+i?3_ + /3+/3_) 



2i?i+i?i_/3+/3_ 



(A.27) 



By substituting Eqs. lfOl) and (|X26l) into Eg. lTOl . 
the optimized fi is given by 



h = 



9T 



80(l + a2)- 



■(2(1 + «^) E E [4*t + 2 (Rl - II)"] + 32a^ E l^l^l + 2 (^?+^t + ^^e)] 

I i=l,3cr=± i=l,3 



(i + 2«2 + «4)E E [$L*L' + 2i?UV] 

CT=±cr'=± 
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64a^ 



(i?l+i?l_ + /i+/i_) (i?3+i?3_ + /3+/3_) + 2i?i+i?i_/3+/3 



(A.28) 



-<5'i+- 



-$2_ 



By minimizing /4 for a fixed ^f_^_+^^_ 
we find that there are two types of quadruple-^ structure, 
i.e., a mixture of four distorted spirals (the ql state) and 
a mixture of two spirals and two SDWs (the q2 state) . 

In the case of a mixture of four spirals (the ql state) , 
/4 takes an extremum, 



(ql) ^ 9T 

^ 80(l + a2)2 



5 + 4a2 + Sa"* 



1 - Sa^ - 14a'' -8a'= + a 



.6 I ^,8 



2(l + a2)' 



(A.29) 



with 



i?i-l 

Ri+ = 
/2 - 



-Rl- — -^3+ 

h- = R3+ 

1 + 4a2 + a^ 
16(l + a2)' ' 
3 + 4a2 + 3q4 



i?3^ 



'+ 16(l + a2)2 



(A.30) 



In this ql state, similarly to the dl state, the four spirals 
are distorted {Ria < ha, ha < Raa) to lower the A2- 
term. For a^ > a^ = (2 - \/3) {ac ^ 0.518), this ql state 
becomes a minimum of fi, while for a < ac it becomes 
locally unstable. 

Next, we consider a mixture of two spirals and two 
SDWs (the q2 state). In the q2 state, a set of #;+ and 
$j_ (z = 1 or 3) form spirals, while the other set of $j+ 
and #j_ (j y^ i) form SDWs perpendicular the spiral 
plane. As an example, we choose here a mixture of two 
spirals characterized by qi± and two SDWs characterized 
by q3± . In this situation, f^ takes an extremum when two 
spirals are isotropic, 

Ri-^h-, (A.31) 

and two SDWs are perpendicular to the spiral plane. 



/3+ 
h- 



(A.32) 



In addition, their amplitudes satisfy 



Ri+ - Ri- 



1 + 12^2 



Ri 



R 



3- 



4(l + 20a2+a4 
4^2 

TTl 

1 + 20^2 + a4 



-m 



(A.33) 



The optimized f^ is then calculated to be 



/. 



(q2) _ 



9T 



40(1+^2)2 



[1 + 12a2 + a*) 
1 + 20a2 + a4 



3 (1 + 4^2 + a*) 



m^. (A.34) 



For a2 < a2 ^ (2 - ^3) [uc ~ 0.518), this q2 state 
becomes a minimum of f^, while for a > ac it becomes 
locally unstable. 

Thus, for the quadruple-g state, we find two metastable 
structures, ql and q2. In FigfTSl we plot /4/r of these 
two metastable quadruple-g structures, given by (JA.29P 
and (jA.34|) . respectively, as a function of a. Note that 
the critical value ac — 0.518 (q;2 = 2 — \/3) of the stability 
is common between the ql state and the q2 state. Since 
a is in the range of a ~ 0.4 for the realistic parameter 
value of J2/J1, the structure characterized by a mixture 
of two spirals and two SDWs (the q2 state) is the stable 
one between the two metastable quadruple-q structures. 



4. The sextuple-g state 

We finally discuss the case where all six wavevectors 
coexist, i.e., the sextuple-q state. In this situation, al- 
though almost all types of combinations of spirals and 
SDWs turn out to be unstable, a mixture of six SDWs 
can be locally stable. 

For a mixture of six SDWs, f^ is minimized when 
SDWs of ql' and q^ run along a common axis while 



% 



S 



SDWs of qf and qf {i ^ j) are mutually orthogonal. 



* 



i± 



*, 



±e 



(1 = 1,2,3). 



(A.35) 



As in the case of the quadruple-g state, the phase factors 
are optimized to minimize the yl4-term. By substituting 



Eq. (IA.35|) into Eq. (IA.5I) , the ^4-term is calculated as where A4 is given by 
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A4 = cos {01+ + 01- - 63+ + 03-) $i+$i_$3+$3_ 

+ COS (6*2- + 02+ - O1+ + 01^) $2+$2-'Sl+'Sl- 
+COS {03^ + 03+ - 02+ + 02-) $3+ $3- "52+ $2-- 



Then, the A^-tevTn is minimized for 

COs{0i+ + 0i--03++03-) = -^, 
COS [02+ + 02- - O1+ + 01-) = -1, 
COS (03+ + ^3.-^2++ 02-) = -1, 



(A.36) 



A. 



$l + $l_$3+$3_ - $2+*2-$l + $l- 



$3+$3_$2+*2-- (A.38) 



From Eqs. (|A.35[) and (jA.37p . the optimized /4 is given 
(A.37) by 



/4 = 



9T 



80(l + a2)2 



12(l + a4)^(<i>4^ + <i>4_)+96«2^<j,2^^2_^8(l^2«2 + a4)^ ^ ^2^ ^2^^ 

i i (ij) cr,cr'— +,— 



r 



(A.39) 



If we fix the total amplitude m^ — ^^ ^^ <i>|^, /4 is min- In this situation, /4 of the sextuple-q state is reduced to 
imized when the amplitudes of six SDWs are in common. 



o m 
$2 ^ 

6 



(A.40) 



,(sc) 

Ji 



3T 



40(1 + a2) 



2^2 



(7 +12^2 + 7a*) 



4\ 4 



(A.41) 
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